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This paper defines an angular velocity for time-dependent functions on the sphere, and applies it to gravitational 
waveforms from compact binaries. Because it is geometrically meaningful and has a clear physical motivation, 
the angular velocity is uniquely useful in helping to solve an important — and largely ignored — problem in models 
of compact binaries: the inverse problem of deducing the physical parameters of a system from the gravitational 
waves alone. It is also used to define the corotating frame of the waveform. When decomposed in this frame, the 
waveform has no rotational dynamics and is therefore as slowly evolving as possible. The resulting simplifications 
lead to straightforward methods for accurately comparing waveforms and constructing hybrids. As formulated in 
this paper, the methods can be applied robustly to both precessing and nonprecessing waveforms, providing a 
clear, comprehensive, and consistent framework for waveform analysis. Explicit implementations of all these 
methods are provided in accompanying computer code. 



I. INTRODUCTION 

Gravitational-wave astronomy stands on the brink of deliv- 
ering numerous observations of merging compact binaries [1- 
6]. Though the uncertainties are large, black-hole binaries 
involving large spins are expected to constitute a significant 
fraction of observable events [7-10]. If these spins are 
misaligned with the orbital angular velocity, the system will 
precess, imprinting the gravitational radiation with strong 
variations [11-14]. While it is not clear how common such 
misalignment will actually be, it is entirely clear that we will 
need good models of the precessing waveforms if we hope to 
accurately measure them. 

We can describe the motion of a precessing binary on 
very short timescales as a simple orbit in a plane; on longer 
timescales, that plane rotates. Now, we know that the 
gravitational-wave field of a nonprecessing system can be 
decomposed into relatively simple modes when the orbital 
plane is orthogonal to the z axis [15-17]. But precession moves 
the orbital plane out of alignment, causing the modes to mix 
and leading to complex behaviors which complicate analysis 
of the waveforms [18-29]. In particular, none of the methods 
developed to analyze nonprecessing systems will work correctly 
with precessing systems. 

In the context of post-Newtonian models, Buonanno, Chen, 
and Vallisneri [18] proposed a convention whereby effects of 
precession can be isolated from orbital motion. Specifically, 
the system is analyzed at each instant in a frame with its z axis 
orthogonal to the orbital plane; from moment to moment, the 
frame is made to rotate to follow the precession. This method 
was later rediscovered in the context of numerical relativity, 
and techniques were developed for finding such a frame from 
the waveform itself in a geometrically meaningful way [30-32]. 

This paper extends previous work by developing a frame 
in which all rotational behavior is eliminated, simplifying the 
waveform as much as possible, and allowing direct generaliza- 
tions of methods for analyzing nonprecessing systems. In the 
process, the angular velocity of a waveform is introduced, which 
also has important uses, such as supplying a partial solution to 
an important inverse problem. 



A. The modeler's inverse problem 

We might distinguish two significant inverse problems 
related to gravitational waves: the modeler's inverse problem 
and the equally important astronomer 's inverse problem. The 
gravitational- wave astronomer's task is to deduce the param- 
eters (masses, spins, etc.) of a system from observations at a 
single point over an extended time. In practice, it is greatly 
complicated by the presence of noise in the data. Usually 
referred to as parameter estimation, this problem has been 
extensively studied [33-44]. The modeler's task, on the other 
hand, is to deduce the parameters given observations of the 
entire sphere at infinity over a brief (possibly infinitesimal) 
interval of time. It is — in some sense — prior to the astronomer's 
problem, because it addresses the meaning of the parameters in 
models astronomers use. 1 This paper concerns itself exclusively 
with the modeler's inverse problem. 

Various methods exist for producing gravitational 
waveforms — numerical-relativity, phenomenological, post- 
Newtonian, and effective-one-body models, for example. 
But no one of these is capable of producing an accurate and 
complete waveform on its own. Numerical-relativity (NR) 
simulations are too expensive to simulate more than a short 
portion of the waveform near merger. Phenomenological 
models use NR data as inputs. Post-Newtonian (PN) 
approximations break down before the merger. Even terms 
generating the effective-one-body (EOB) inspiral and the ad 
hoc method of attaching a ringdown must be "calibrated" by 
comparison to numerical results. Therefore, we need more 
than one model to generate a complete waveform, which means 
that we need to understand precisely how the different models 
relate to each other. This leads directly to the inverse problem. 

The numbers we plug in to a computation of initial data for 
an NR simulation bear no clear relation to the numbers we 
insert into a PN or an EOB computation — or even to other 
NR simulations using different formulations. For example, 



Intriguingly, understanding the modeler's inverse problem may help to 
inform the astronomer's inverse problem more directly [29]. 
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the direction of a black-hole spin measured in the arbitrary 
coordinates of NR initial data need not correspond in any 
meaningful way to the direction measured in PN coordinates. 
Even if the gauge condition used for the numerical simulation 
were the same as the one used to derive the PN formulas, the 
initial data would not be the same, so the gauge itself would 
be different. For nonprecessing systems, symmetries reduce 
the ambiguity to one of simple time and phase offsets, and 
numerous simple methods have been suggested to resolve those 
ambiguities [45]. But in the precessing case, we need to be 
much more careful. Simply using the same numbers in two 
different models leads to comparing systems with inherently 
different physics. 

Fundamentally, we need to establish a mapping between the 
input parameters of different models such that they produce the 
same physics (as nearly as possible) during some span of time 
for which both models are valid. Because we have no access to 
any invariant physical meaning behind our parameters, we need 
to take a different approach. For example, given some particular 
set of parameters, we can run a numerical simulation. Then, 
we can work backwards from the resulting waveform and try 
to find the parameters needed to generate the same waveform 
with a PN system — which is the inverse problem. 

The issue of ascribing meaningful physical interpretations 
to geometric quantities measured on has been investigated 
to some extent [46-48], but it is not clear that these methods 
are useful for the immediate problem of analyzing gravitational 
waveforms. The angular- velocity vector introduced by this 
paper and a related vector introduced by O'Shaughnessy 
et al. [31] provide geometrically meaningful physical quantities 
which can be measured directly from the waveforms alone, and 
are thus prime candidates for use in solving the inverse problem. 
Indeed, we will see in Sec. Ill that these two vectors are very 
closely related to input parameters for the precessing PN system. 
This provides a partial solution to the inverse problem, leaving 
three remaining degrees of freedom. Several possibilities will 
also be suggested for completing the solution of the inverse 
problem, though they are beyond the scope of this paper. 

B. Overview of this paper 

Section 11 introduces the angular velocity to of a waveform, 
finding a straightforward formula and a more intuitive inter- 
pretation of the mathematics. This and the related vector 
Vf suggested in Ref. [31] are then used in Sec. Ill to find a 
partial solution to the inverse problem. It is shown that co and 
V f — which are measured from the waveform alone — can be 
combined to give expressions for the corresponding PN orbital 
elements. This can then be used to deduce the parameters of 
the system. A PN waveform is used as a test case, showing 
excellent agreement between the original parameters and the 
parameters deduced from the waveform alone. In Sec. IV, 
the angular velocity is used to determine a frame with that 
velocity. The same PN waveform used in the previous section 
is decomposed in this frame, showing that the amplitudes of the 
waveform modes become very simple, and their phases become 



nearly constant. Because this frame reduces the complexity of 
the waveform, it is ideally suited to practical manipulation of 
waveforms, which is discussed in detail in Sec. IV C. It is also 
worth noting that the partial solution to the inverse problem 
completely establishes all extrinsic parameters, giving us a 
solid foundation for comparisons between waveforms. Finally, 
the results are summarized and suggestions for future work are 
collected in Sec. V. 

The appendices provide deeper background information 
which may be useful for implementing these methods or 
comparing to other methods. Appendix A presents a fairly 
comprehensive discussion of quaternions and various related 
details, including several new results. In Appendix B, formulas 
are derived for the rotation of arbitrary spin-weighted functions. 
While equivalent formulas have been derived previously [32, 
49, 50], this derivation uses a somewhat different technique, and 
carefully develops conventions for consistency throughout this 
paper. In any case, the upshot is that modes of spin-weighted 
fields transform exactly as do modes of spin-weight-zero fields. 
Lastly, Appendix C discusses related previous work in the 
same formalism used in this paper, allowing for more direct 
comparisons. 

Ancillary files included with this paper (available on the 
paper's arXiv page) contain computer code implementing all of 
the concepts introduced here, among others. The core functions 
are written in C++ [51] for speed, using several functions from 
the GNU Scientific Library [52]. While this code could be 
incorporated directly into other C/C++ codes, an additional user 
interface is provided in Python [53, 54] code as the GWFrames 
module, which simply exposes all the C++ functions through 
Python. Documentation and examples can be found among 
the ancillary files. Relevant functions or classes are mentioned 
where appropriate throughout this paper. 

C. Quaternion notation 

The techniques of this paper necessarily involve rotations, 
which are best implemented in terms of the group of unit 
quaternions because of the numerous advantages over direct 
manipulation of rotation matrices or Euler-angle coordinates. 
By using quaternions, we obtain robust methods that can 
be blindly applied to general systems, including nonprecess- 
ing ones — which simplifies the processing of large numbers 
of waveforms. Moreover, once the basics are understood, 
quaternion rotations are more intuitive than either of those 
inferior descriptions. In fact, quaternions are essentially the 
axis-angle description of rotations, in a more practical guise. 
Therefore, quaternion notation will be used throughout. As 
mentioned above, Appendix A provides a thorough introduction 
to quaternions, while computer code included in ancillary files 
with this paper gives practical implementations of the necessary 
functionality through GWFrames. Quaternion. However, such 
details are not necessary for a good understanding of this paper; 
the following paragraph should provide sufficient background. 

Quaternions can be thought of as generalizing the familiar 
complex numbers, where the imaginary part is generalized to a 
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three-dimensional vector part. 2 We can write a quaternion as 
the sum of a scalar and a vector: Q = qo + q. The conjugate 
of the quaternion is Q = qo - q. We can multiply quaternions 
together [see Eq. (A2)], the product being associative but not 
commutative in general. The norm of a quaternion is defined 
according to |Q| 2 = QQ = + q • q. Unit quaternions, 
having norm |R| — 1, are especially important, as they describe 
rotations. To see this, we can consider a vector to be a 
quaternion with scalar component equal to zero, in which case 
it makes sense to multiply a vector by a quaternion. Then we 
can define the transformation 



RvR 



(1) 



A simple exercise shows that this transformation is linear, 
and preserves lengths and orientations, so it is just a rotation. 
Ultimately, the best reason to use quaternions is the existence of 
simple formulas [Eqs. (A8) and (A9)] for the exponential and 
logarithm, which prove to be endlessly useful. In particular, 
we can express an arbitrary unit quaternion as R = e 6 "/ 2 = 
cos | +u sin I, where exponentiation of a quaternion is defined 
by the usual power series. 1 It turns out that this R produces 
a rotation through the angle 9 about the axis u. Because 
any rotation may be expressed in this form, we will use unit 
quaternions as our only representation of rotations, and refer 
to them as rotors. Conversely, given a rotor R, we can find the 
corresponding axis and angle according to 6u =2 log R, where 
the logarithm is given by Eq. (A9). A frame will be described 
by the rotor that generates it by rotating some standard basis 
frame. 



II. ANGULAR VELOCITY OF A WAVEFORM 

We can define the angular velocity of a gravitational 
waveform — or any field on a sphere — as the opposite of the 
velocity of the counter-rotation needed to keep the field as 
constant as possible. In the first part of this section, this 
definition will be formulated more precisely, resulting in a 
surprisingly simple formula for the angular velocity. The 
formula can be interpreted as a projection of the familiar 
operator equation -ico ■ L - d t onto the "rotational parts" 
of the waveform — a notion which can be made surprisingly 
rigorous using the language of Hilbert spaces, as discussed in 
the second part of this section. 

A. Finding the angular velocity 

The essential idea here is to remove the rotational behavior 
of the waveform by imposing a rotation that eliminates as much 



2 In fact, both complex numbers and quaternions are special cases of geometric 
algebra [55], done in two and three dimensions, respectively. Much of 
our intuition from complex algebra transfers easily to quaternion algebra 
when i = V— T is replaced by a unit vector. The notable exception to this 
correspondence is noncommutativity of the quaternion product. 

3 Note the striking — and not coincidental — similarity to Euler's formula with 
u in place of the unit imaginary i. This results from the fact that, under 
quaternion multiplication, «« = -!. 



of the time dependence as possible. Suppose that Mj(t) is a 
time-dependent rotation operator acting on the wave field / 
(usually representing ¥4 or h) such that Mj(tj) = 1 . We wish 
to find the rotation operator that — in some sense — minimizes 
the quantity 



dt 



[#,(*) /ft ft?)] 



(2) 



Clearly, this is a complex function of position on the sphere. To 
reduce it to a single real number, we take its squared magnitude 
and integrate over the sphere: 



S(^):= ^ |[^(f)/(f;#,0] 



dQ 



(3) 



If we expand the field / in spin-weighted spherical harmonics 
(SWSHs, discussed in Appendix B), the natural way to express 
the rotation operator is in its usual form 3$j = exp[-i Oj ■ L], 
where L is the standard angular-momentum operator and Oj 
is the time-dependent axis-angle description of the rotation. 
Note that we must have Oj(tj) = because we have assumed 
that &j(tj) = 1 . This is absolutely crucial because it makes the 
differentiation in Eq. (3) tractable. We also define the angular 
velocity 4 



co := -d,8j 



(4) 



where the negative sign arises because 0j corresponds to the 
rotation needed to keep the field fixed in a moving frame, 
whereas co is intended to describe the motion of the field relative 
to the initial static frame. We have 



B(<»)= f 

Js 2 



\ico- Lf + d,f\ 2 dO. 



(5) 



Now, we can write the integral in terms of a sum over standard 
matrix elements of the angular-momentum operator and the 
problem simplifies nicely. We obtain 



S = co ■ (LL) ■co + 2co-{L8 t } + Y J \d t f 



e,m 



where we have defined the matrix 5 

{LL) ab := J] f' m ' {{,m'\L (a L b) \{,m) f' m , 



(6) 



(7a) 



4 Subscripts are necessary on the rotation operator Mi and the associated 
vector 6j because these have certain properties depending on which instant 
of time tj we are looking at. We have implicitly assumed that the Oi(i) are 
all related by simple constant offsets, as necessary to satisfy the conditions 
Oj(tj) = 0. Because the offsets are constant, the angular-velocity vector oj 
does not have such a dependence, and so does not need the subscript. 

5 Here, the \l, m) represent the spin-weighted eigenfunctions, but the angular- 
momentum operator acts on these just as in the non-spin- weighted case [49], 
making this notation particularly familiar. The matrix denoted here as {LL) ab 
is precisely the quantity (L( a Z,m)i defined by O'Shaughnessy et al. [31], 
except that the latter is normalized by 2f,ml/'' m | 2 ' 
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and the vector 

{Ld,) a := J] 3 [f' m ' {(,m'\L a \{,m) d,f' m \ . (7b) 

£,m,m' 

Noting that the last term in Eq. (6) is independent of to, we can 
find the minimum 6 of S analytically: 



to = -{LLy x ■ {Ld,) 



(7c) 



The effects of L a are familiar, so this may be directly computed 
from knowledge of f e ' m (t), with no optimization or solution 
of the eigensystem necessary. There is no ambiguity in the 
direction of the angular velocity, and we obtain a meaningful 
magnitude. 

In the computer code included among this paper's an- 
cillary files, a waveform object may be constructed with 
GWFrames. Waveform. The angular velocity may then be found 
using the AngularVelocityVector method on such an object. 

B. Interpreting the mathematics 

Equation (7) gives a formula for the angular-velocity vector 
of the field /. Though it takes a relatively simple form, the 
reason it takes this particular form may seem somewhat opaque. 
In fact, it really has quite a simple interpretation, which may be 
instructive. In fact, we can start off with a simple observation 
and re-derive Eq. (7) in a very different way. 

The angular-momentum operator L generates rotations, as 
is well known. So, for example, -\co ■ L f gives the time 
rate of change of the field under a simple rotation given by 
to. More generally, the three components of -i L f form a 
basis generating the Hilbert subspace A, consisting of functions 
describing possible rates of change for / under (complex) 
rotations. On the other hand, we also have a second operator d,, 
which gives the actual time rate of change of the field, whether 
that change is a simple rotation or a change in amplitude — or 
a more complicated behavior. But we can extract the part of 
d t f caused by (real) rotation alone by projecting onto the basis 
vectors of A and taking the real part: 



ilr- 



iLfdJdQ 



(8) 



The three components of this expression completely describe 
the rotational part of d t f. We take the real part because we 



' We can show that it is a true minimum rather than a more general stationary 
point by looking at the Hessian matrix of 5, which is just 2 (LL). We are free 
to rotate this matrix into a frame in which its dominant principal axis is along 
Z. Then, we can calculate its eigenvalues as (L 2 ) and (L 2 — L 2 ) ± \(L 2 )\. 
As long as some mode with m # is nonzero, these are always (strictly) 
positive. Hence, (LL) is positive definite, and we have a true minimum. 
Furthermore, we can calculate that the determinant is actually the product 
of these eigenvalues, and thus is also nonzero whenever the field is nonzero, 
allowing us to invert the matrix in Eq. (7c). Since (LL) is a geometric 
object, and eigenvalues and determinants are invariant under rotations, these 
conclusions hold in all frames. 



ordinarily take the dot product of -i L with a real-valued vector, 
so if we expect to find such terms in d t f, they must have real 
components [56]. 

Now, the crucial point: if to correctly describes the rotation, 
the same projection of -i to ■ Lf must give the same result: 



-\Lf(-\to-Lf) dQ. 



-iL/<9,/dQ 



(9) 

If we expand / in spin-weighted spherical harmonics, it turns 
out 7 that this equation reduces to precisely 



{LL)-to = -{Ld t ) , 



(10) 



which is, of course, equivalent to Eq. (7c). 

Thus, we see the interpretation clearly. In the case of a pure 
rotation, we have -i to ■ L f = d,f. In general, however, we 
have to project onto the rotational parts of the waveform for 
equality to hold, which is just what Eqs. (7) and (10) do. It 
is also worth noting that in the purely rotational case, we can 
use -i to ■ L = d, directly and calculate SsO. Recalling the 
definition of S in Eq. (3), this says that the time variation is 
completely eliminated. 

Interestingly, we can see this projection working directly by 
showing that (LL) and (L d t ) are insensitive to changes in the 
amplitudes of the modes. Clearly, (LL) does not depend on any 
derivatives with respect to time. To see that (L d t ) is insensitive 
to changing amplitude, we first note that it is a geometric object 
so we can evaluate it in any frame we choose — we choose a 
frame in which it is aligned with the z axis. Next, we decompose 
the field into (logarithmic) amplitude and phase parts: 



/ f - m (0-exp[/- m (f) + i/'" ! (0] 



(11) 



A pure rotation about the z axis leads to % e,m = and <p e,m = 
-m \co\, so we expect that a projection onto the rotational part 
will eliminate x C ' m but must not eliminate <f> C ' m . In fact, we can 
explicitly calculate 



(L d,) = z J] 3 [/•'" <f, m\L : % m) 8,f 

= 2 2 3 [(**» + i ***) 
e,m 



(12a) 



(12b) 



(12c) 



Here, taking the imaginary part has caused x to drop out 
entirely, leaving only <f> ,m , supporting the claim that we have 
removed non-rotational parts of the waveform. Note that this 
formula is entirely general; we have not assumed any particular 
behavior of^f or 0, for example. 



7 As usual, we get integrals of the form J ...\&,ip) (&, ip\ ... dfi, which are just 
resolutions of the identity. Then, taking the real part on the left-hand side 
is equivalent to symmetrizing over the indices of the two L vectors before 
contracting with &>. On the right-hand side, taking the real part of i times a 
quantity is the same as taking the negative imaginary part, so this is precisely 
the definition of - (L d t ). 
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III. SOLVING THE INVERSE PROBLEM 

Section 1 A established the need to solve the inverse problem. 
Essentially, in order to create a complete gravitational wave- 
form, we need to be able to take a finite or even infinitesimal 
portion of a waveform and infer the parameters of the PN (or 
similar) system that result in that waveform. Because it is 
the most extensively developed system, we will discuss the 
quasicircular PN model as a concrete example. In this section, 
we will first describe the parameters that need to be established. 
This will involve reviewing the basic elements of the PN model. 
We will then see how to solve part of the inverse problem using 
the angular-velocity vector a> and the dominant eigenvector 
of (LL) (denoted V/) [31], showing the effectiveness of this 
method with an example. 

A. The required parameters 

To the extent that different formulations of the PN model are 
correct, they predict the same physics, and so we are free to 
choose between them as we wish. Certain formulations may be 
better than others with regard to solving the inverse problem. 
Here, we follow Refs. [24, 57]. First, we assume a pair of 
particles with masses M\ and M2, and spins Si and 52. The 
unit vector pointing from the second to the first is h. The orbital 
angular velocity is defined as 

il or b := nxh . (13a) 

The direction of this vector is frequently expressed in the 
literature as L N (= O or b). There is (in general) an additional 
rotation of the system due to precession, denoted Q prec . Now, 
if this vector were to have any component orthogonal to h, that 
would contradict the definition of £l or b, so it must simply be 
proportional to n: s 

^prec ■= Oprec « ■ (13b) 

We also define the sum of these: 

fitot := fiorb + fiprec ■ (13c) 

During the evolution we must record the minimal-rotation 
frame 9 aligned with fi or b and the accumulated orbital phase 
O or b measured relative to h. Then, the frame of the binary will 
given by rotating the minimal-rotation frame by (£> 01 -b about its 
z axis. These are the orbital elements of the system. Their 
evolution is not of particular concern here, as the details have 
no effect on our conclusions. The waveform can be calculated 
in this frame using standard formulas, and transformed to an 
inertial frame if needed to complete the construction of the 
waveform. 



Note that for other formulations of the PN model, this equation may not be 
true. See Refs. [11, 12, 18-25] for more details. 
9 See Sec. C4andRef. [32]. 



The initial data we need to begin a PN calculation, then, 
are the values of (Mi , M2, Si, S2, Qorb, n) at some initial time. 
These might be termed the intrinsic parameters of the sys- 
tem [18, 58, 59]. They are geometrically meaningful, and 
covariant under certain symmetries assumed for our system — 
namely time translation and rotation of coordinates. But this 
brings up a subtlety. We can think of two more classes of 
parameters: the extrinsic and the fiducial. Extrinsic parameters 
depend on the observer, and can be thought of generally in 
terms of degrees of gauge freedom like the time offset or an 
overall rotation. Fiducial parameters are selected values of 
intrinsic quantities that depend on extrinsic parameters. By 
solving for the intrinsic parameters relative to a particular time 
function and a particular basis for the vectors, we will be tacitly 
setting the extrinsic parameters. Then, when comparing two 
waveforms, we must choose fiducial parameters and ensure that 
the extrinsic parameters are the same. We will find that it is a 
simple matter to ascertain the intrinsic parameters except for 
three degrees of freedom in the directions of the spin vectors. It 
will also be straightforward to completely establish the extrinsic 
parameters. 

B. Deducing the parameters 

We might only expect quantities to be meaningful if they 
are covariant objects measured at infinity — e.g., waveforms 
or ADM-type quantities. Coordinate locations of black holes 
in a simulation, for example, depend too much on details of 
gauge conditions and vagaries of initial data and junk radiation 
to be of any real use. On the other hand, some quantities are 
also reasonably well defined when the black holes are very 
widely separated. Therefore, if we find that certain quantities 
change slowly during the early part of the NR simulation, 
and are not expected to have changed much previously, then 
we might also be able to use those quantities in our analysis. 
This is typically true of the masses and spin magnitudes when 
measured appropriately [60-62], except to the extent that they 
are expected to change [63, 64]. Therefore, we assume that Mi, 
Mi, |Si|, and can be measured in the simulation and used 
directly. The rest of our intrinsic parameters will come from 
the waveforms (or possibly other measurements on ,# + ). 

To see how we can derive orbital elements from quantities 
observable from the waveform, we need to see how orbital ele- 
ments give rise to the waveform. Familiar calculations [15, 16] 
tell us that the PN waveform is created by motion of the binary. 
The complete motion is described by Q tot , so we expect that oj 
should be the same. On the other hand, the component along 
h does not lead to changing multipole moments (to our level 
of approximation). So only the component of 12 tot orthogonal 
to h is involved — but that is precisely £l OI b- We can therefore 
expect that the waveform is oriented along this vector, in some 
sense. Now, the vector z happens to be the dominant eigenvector 
of (LL) for an individual spin-weighted spherical harmonic 
(though not necessarily for a combination of them). It also 
happens to be the dominant eigenvector when the field / is 
symmetric under reflection through the x-y plane [65] — as the 



5 



MICHAEL BOYLE 



PN waveform is in the frame aligned with S2 rt>- Therefore, we 
should expect fi or b to be parallel to Vf. 

Putting these considerations together, we can expect the 
following approximate equalities: 



<jj . 



n orb - (Vf -io)v f , 
il prec =i h) - \Vf ■ a>) V 



f 



(14a) 
(14b) 
(14c) 



Inspection of the PN model suggests that these expressions 
should become more exactly true in the asymptotic limit of low 
orbital velocities. Note that Eq. (13b) shows that h is along 
i^prec, so we effectively obtain that quantity as well, whenever 
the precession is nonzero. 

Figure 1 compares the orbital elements to the related wave- 
form expressions, for a PN system with significant precession. 
The direction of £2 or b coincides extremely well with Vf — they 
agree to within the numerical precision throughout the inspiral. 
fitot an d w agree to within a few parts in 10 5 early in the inspiral, 
though the disagreement grows near merger. However, it may 
be possible to remove even this disagreement through more 
careful treatment of the distinction between the orbital phase 
and the phase of a waveform in PN theory. 

Now, in each case of Eqs. (14), the quantities on the right- 
hand side are measured directly from the waveform. Thus, 
if we have a numerical waveform, we can simply measure 
the right-hand sides and define the "PN-equivalent" orbital 
elements according to these equations. A PN system given 
those parameters as initial conditions will necessarily be as 
similar to the numerical system as possible — at least by the 
measures of a> and Vf. 

This gives us a partial solution to the inverse problem. We 
are lacking four degrees of freedom corresponding to the 
directions of the two spin vectors. One additional piece of 
information is also available from the foregoing considerations. 
The magnitude Q prec is given in PN theory as a bilinear function 
of Si ■ h and 52 • h, where the coefficients depend on the 
PN-expansion parameter v := (MQ or b)^ 3 , and are therefore 
already known. Thus, we can solve for Si • ft, for example. 
This suggests other possible methods to find the remaining 
three components of spin. For example, PN expressions for 
the orbital angular momentum L are available in terms of the 
orbital elements and the various projections of the spin vectors. 
If it is practical to measure the total angular momentum of the 
spacetime J in the numerical solution [66, 67], we could then 
use the PN expression for L + S to solve for the PN-equivalent 
components of spin. This would complete the solution of the 
inverse problem. 

Alternatively, we might measure various modes of the wave- 
form and equate them to the PN expressions for those modes. 
Again, these expressions contain various known quantities, as 
well as bilinear combinations of Si ■ £2 01 -b and S2 • £l OT b [13, 68]. 
Therefore, we could solve for these combinations of the spin 
components. Seemingly, this would rely on the accuracy of 
the PN expressions, which is not very high for spin terms. On 
the other hand, the influence of any errors that result would 



be similarly diminished. One final degree of freedom would 
remain in this example, and would have to be fixed by other 
means. In any case, we leave these considerations to future 
work. 

In the computer code included among this paper's an- 
cillary files, a waveform object may be constructed with 
GWFrames. Wave form, or a PN waveform may be constructed 
with GWFrames .PNWaveform. The PN-equivalent orbital and 
precessional angular velocities may then be calculated using the 
PNEquivalentOrbitalAV and PNEquivalentPrecessionalAV 
methods. 

IV. THE COROTATING FRAME 

So far, we have calculated only the instantaneous angular 
velocity of the waveform, a> [Eq. (7c)]. While this has 
already proven useful in the previous section, it can also be 
advantageous in determining a frame in which to decompose the 
waveform. Specifically, we seek a frame whose angular velocity 
is just a>. When decomposed in this frame, the waveform will 
have no rotation, and will be as constant as possible. The frame 
compares favorably to other frames introduced previously [30- 
32] (see Appendix C). It has practical benefits and suggests 
simple techniques for measuring, comparing, and processing 
waveforms from numerical simulations. 

A. Finding the corotating frame 

Our task here is to find the rotor R(f) describing a frame 
whose angular velocity is a>. We can relate the two by a simple 
equation: 



w=2RR. 



(15) 



(See Ref. [32] and Sec. A3.) Unfortunately, the solution we 
might naively write down is wrong: 



R(f) + exp 



a>(t') dt' 



(16) 



except when the system is nonprecessing. Ultimately, the reason 
for the failure of this formula in general is that co is not parallel 
to its derivative (or integral). In the language of quaternions, 10 
a> fails to commute with its derivative (or integral); to find 
the correct version of Eq. ( 1 6), we need to account for that 
noncommutativity. 

Given a>, we could solve Eq. (15) for R and integrate as we 
would a vector equation. But in practice this would quickly 
violate the constraint that R should be a unit quaternion. Instead, 
we will need an expression in terms of the logarithm of this 
rotor: 



RR 



r + 



sin 2 |t| 
2|c| 2 



M]+ W-mWco, W 
4 r 3 



' Note that the failure to commute is by no means specific to the quaternion 
description of rotations; it is a feature of rotations themselves. Quaternions 
do, however, provide a very effective means of solving the problem. 
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FIG. 1. Orbital elements compared to waveform quantities. These plots show the PN orbital elements £l olb (left) and Sl lol (right), compared 
to the "PN-equivalent" quantities derived from the waveforms alone given in Eqs. (14). The binary has a 6 : 1 mass ratio. Initially, the larger black 
hole has a spin of 5 i/M 2 { = 0.9 in the (i?, <p) = (2.00, 0.25) direction; the smaller black hole has a spin of S 2 IM\ = 0.3 in the (■&, if) = (2.4, 2.9) 
direction. These parameters were chosen because the resulting orbital velocity happens to execute a complete flip, passing very close to —z, 
which is a rigorous test of these methods. In both plots, tighter oscillations correspond to earlier times; the last 3600 M before merger are shown. 
The directions Ci mb and Vi, are identical to within the numerical accuracy throughout the inspiral. The vectors £l tot and a> are the same to within 
a few parts in 10 5 early in the inspiral, though differences grow somewhat as the system approaches merger (roughly the end of the data shown 
here). 



with r(f) := log R(f). (See Sec. A 3 for the derivation.) The 
second and third terms on the right-hand side of this expression 
account for the noncommutativity as needed. Setting the right- 
hand side equal to cj/2, we can solve for r to find 



t(t-o})\ Irl cotlrl r(r • cS) 1 

t = \(D „ — — + — ^ + -CJXX. 

2 2 Irl 2 2 



(18a) 



This is just an explicit first-order ordinary differential equation, 
so we can integrate numerically using standard techniques to 
arrive at the appropriate r(f) and find the corotating frame 



RcorotW = exp [r(f)] 



(18b) 



Using the fact that R = exp r and -R = exp V^gr *] describe 
the same frame, we can reset the value of r between steps of 
the numerical integration to keep its magnitude small. This 
improves the quality of the numerical integration, though it 
may then be useful to go back and flip the signs of rotors as 
necessary to keep R(?) as continuous as possible. The procedure 
is described in more detail in Appendix A3. 

The advantage of this method over direct integration of 
Eq. (15) is that it ensures that the resulting quaternion truly 
does have norm 1 . When integrated directly, the quaternion 
in Eq. (15) has four degrees of freedom, whereas a unit 
quaternion has only three. By integrating Eq. (18a) instead, we 
eliminate the extra degree of freedom, reducing this to a truly 
three-dimensional problem while automatically satisfying the 
constraint on the norm. In general, transforming equations 
in such a way improves the accuracy of numerical results 
significantly — as is certainly the case with this system when 
tested. 



Naturally, imposing a condition on the angular velocity of 
a frame leaves its overall orientation free. Assuming R CO rot(0 
describes a frame whose angular velocity is a>, then the frame 
R CO rot(0 Rc will hare the same angular velocity for any constant 
R c . Alternatively, the frame R c R CO rot(0 would have angular 
velocity co rotated by R c . In the interests of simplifying the 
waveform, it is best to choose some particular time during the 
inspiral at which to align the z axis of the frame with V/, as 
suggested by O'Shaughnessy et al. [31]. Once this is done at 
one instant of time, z and Vf should be aligned at all other 
times to very high accuracy. We are still free to rotate about the 
z axis, so we can set the phase of the ({, m) = (2, 2) mode to 
at this instant, for example. We will see below that the phase is 
very slowly varying in the corotating frame, so this will not be a 
delicate operation. Alternatively, if the waveform is precessing, 
we can align the x axis with the PN-equivalent h, which should 
be roughly equivalent to setting the (2, 2) phase to 0. When 
comparing two waveforms, the only requirement is that these 
instants of time be comparable, which is assured by choosing a 
common fiducial quantity. These issues are discussed further 
in Sec. IV C. 



In the computer code included among this paper's ancil- 
lary files, the function GWFrames.FrameFromAngularVelocity 
returns the corotating frame, given an array of quaternions 
representing the angular-velocity vector as a function of time. 
Alternatively, a waveform object may be constructed with 
GWFrames. Wave form and transformed to the corotating frame 
with the TransformToCorotatingFrame method. 
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B. Gravitational waveforms in the corotating frame 

Figure 2 demonstrates the effects of decomposing the PN 
waveform described in Fig. 1 in various frames. First is the 
usual inertial frame, where a stationary observer at infinity has 
constant coordinate position. In this frame, the moduli of the 
modes oscillate wildly, as power shifts between them (upper 
left panel of Fig. 2). Similarly, the phase (upper right panel) 
shows strange features. The I = ±2 and I = ±1 modes have 
roughly the same frequency, as power from the dominant modes 
leaks into and overwhelms the € = ±1 modes. Those phases 
also change direction each time the rotation axis passes through 
the x-y plane (at times of roughly -2400 and -600). Both 
the modulus and phase are very complicated functions in this 
frame. They would be hard to model directly, and their rapid 
variations are not conducive to accurate numerics. 

The next pair of panels shows the waveform in the waveform- 
aligned frame suggested by O'Shaughnessy et al. [31] supple- 
mented with the minimal-rotation condition (see Sec. C 4 and 
Ref. [32]). Here, (LL) is evaluated using all modes up to t — 8, 
the dominant eigenvector is found, and the fame is rotated so 
that its z axis coincides with that eigenvector while obeying the 
minimal-rotation condition. This drastically simplifies both the 
modulus and phase, as seen in the middle panels of Fig. 2. The 
modulus is very smoothly sweeping up as the binary spirals in 
toward merger. The phases are now separated as usual, with 
slopes more nearly equal to — nO^. 

Decomposing the waveform in the waveform-aligned frame 
also requires recording the orientation of the frame. In that case, 
there is no additional overhead in going to the corotating frame, 
which is shown in the lower panels of Fig. 2. The modulus plot 
is identical to the one in the previous case. However, there is 
further improvement in the phase, with each mode having nearly 
constant phase throughout. Similar behavior is seen in other 
modes with t > 2. Naturally, such a waveform is particularly 
well suited to interpolation and hybridization [26, 32, 72]. 

A minor feature to note in the phase is the non-constancy of 
the (2, +1) modes. Considered on their own, these variations 
could be removed by a rotation because the curves change in 
opposite directions. On the other hand, this would cause the 
phases of other modes to vary. Equation (7c) automatically 
balances these concerns; the amplitudes of the (2, +1) modes 
are so small that they do not carry much weight. By transform- 
ing to the corotating frame, we isolate the waveform's intrinsic 
dynamics — seen here in 2 ±1 — from the rotational dynamics 
of the system, allowing for separate analyses. This is important 
because they are separately modeled, so it is useful to be able 
to inspect each effect on its own. 



C. Extrapolation, comparison, alignment, and 
hybridization 

As a practical matter, we need to manipulate numerical 
waveforms in various ways. We must eliminate physical and 
gauge effects associated with extraction of data at finite radius, 



usually by extrapolation [73-75]." To compare numerical 
waveforms to each other or to analytical waveforms, we need to 
determine the extrinsic parameters corresponding to freedom in 
choosing the zero of time and the overall orientation of our axes. 
To construct a complete waveform, we may occasionally need to 
hybridize waveforms from different systems [26, 32, 72, 79-82]. 
Various approaches to these problems have been introduced, 
but most are designed exclusively for nonprecessing systems, 
or are otherwise incomplete or fragile. The angular velocity 
and corotating frame provide excellent tools for addressing 
these issues more generally, and are especially robust when 
implemented by means of quaternions. 

Practical extrapolation relies on smoothness of the extrapo- 
lated functions [75]. In the corotating frame, the modes of the 
waveform are essentially constant during inspiral — except for 
the overall growth in modulus — suggesting that this is the ideal 
frame for extrapolation. We can implement such a procedure in 
the usual way, with one minor addition. We first impose a time- 
retardation offset to all the data, as usual. Then, we add the 
step of finding the corotating frame of the outermost extracted 
data, and transforming the data at all radii to that frame. (We 
cannot rotate data at each radius into its own corotating frame, 
as that would require extrapolation of rotors, which is not well 
understood.) The waveforms at the other radii will not be 
precisely in their own corotating frames, but should be close 
enough that the data are quite smooth. The extrapolation may 
then proceed as usual, resulting in an extrapolated waveform 
in the corotating frame of the outermost data. Again, the 
extrapolated result will not be precisely its own corotating 
frame, but the transformation is routine. 

In the nonprecessing case, a multitude of methods have 
been suggested to compare and align waveforms (fixing the 
extrinsic parameters) and to construct hybrids of NR and PN 
waveforms [45]. These all require generalization to use in 
the case of significant precession. By using the corotating 
frame, we simplify the modes of the waveform decomposition 
sufficiently that ordinary methods can still be used. (Com- 
paring phase differences or relative differences in modulus, 
for example.) But each waveform now comes with its own 
frame, Ra(0 and Rfl(f), encoding most of the phase dynamics, 
so we will also need to compare the difference between the 
frames themselves. Fortunately, quaternions provide us with 
geometrically meaningful measures of the difference. 

The difference itself is given in quaternion form as 12 

R A (t):=R A (t)R B (t). (19a) 



1 1 Cauchy-characteristic extraction is another method of finding the correct 
waveform at [76-78]. This can proceed as usual, and the final waveform 
can be transformed to the corotating frame. A recent implementation [78] 
has improved the efficiency to make this a more attractive alternative to 
extrapolation. 

12 The inverse of an arbitrary nonzero quaternion Q is just Q/|Q| 2 . Since rotors 
have norm 1 , the inverse of a rotor R is just R. Therefore, this formula is 
analogous to subtraction, but applied to rotation operators. 
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FIG. 2. Precessing waveform in various frames. These plots show the modulus (left) and phase (right) of the I = 2 modes of a post-Newtonian 
waveform in the inertial (top), waveform-aligned minimal-rotation (middle), and corotating (bottom) frames. The phase is defined as usual [69- 
71] so that //'" = |//-"'| exp[i0 f "'], with branch-cut discontinuities removed. The system is the same as the one shown in Fig. I. Going from the 
inertial frame to the aligned frame drastically simplifies the waveform amplitudes and significantly simplifies the phase. In the aligned frame, the 
waveform looks very much like a nonprecessing waveform [26]. Expressing the waveform in the corotating frame retains the smoothness in 
amplitude seen in the aligned frame, but makes the phases of the modes nearly constant, with values of roughly 0, ±n/2, and n. Similar results 
can also be seen for modes with t > 2. 



This is the rotation taking frame B into frame A, and is 
independent of the basis with respect to which A and B are 
defined. Now, we might want to know how "big" this difference 
rotation is. As noted in Sec. IC, we can write any rotation, 
including Ra, in axis-angle form: 

R A (0 = exp[O A /2] . (19b) 

We can easily solve this equation for Oa by taking the logarithm. 
In particular, its magnitude is the angle through which the 
system must be rotated: 

<D A (*)=2|logR A (*)| . (19c) 

This can be used as a simple but complete description of the 



phase difference between two systems. 

We can understand this better and make contact with previous 
work by recognizing Oa as a more general version of a common 
measure of the difference between waveforms common in 
analysis of nonprecessing systems. That measure is A(f> t,m , the 
difference between the phases of the modes as measured in the 
static frame. Because the angular velocity is conventionally 
chosen to be along the z axis, we can usually relate the orbital 
phase <5 01 -b to the waveform phase as <p l ' m as -m <D or b. Here, we 



It is crucial to note that log(R,i Rg) # logR.4 - logRg because rotations 
do not commute. In particular, the latter depends on the basis frame, and is 
therefore not a useful measure of the difference between frames. 
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have the similar expression 

A0°" « -m<E> A . (20) 

However, as we saw in the previous section, A<p e,m is a less 
useful measure for precessing systems. In contrast, Oa encap- 
sulates the differences in both the orbital and the precessional 
dynamics 14 in one convenient function while leaving the 
waveform dynamics separate, and is equally relevant in both 
precessing and nonprecessing systems. 

The quantity €>a gives us a compact description of the 
difference between two waveforms as measured in their coro- 
tating frames. But it depends on the extrinsic parameters 
discussed in Sec. Ill A: the overall time offset and orientation 
of the static basis frame. Alignment of waveforms consists 
of minimizing differences between the waveforms at some 
instant or over some span of time by adjusting the extrinsic 
parameters as needed. This can be seen as a restricted 
version of the inverse problem, where we simply assume that 
the intrinsic parameters are identical — as when we wish to 
compare waveforms evolved the same initial data with different 
numerical resolution. Section 111 B used implicit assumptions 
that the time coordinate and basis frame would be defined to 
be the same in both waveforms, and the physical parameters 
are expected to be the same in that frame. Here, we are simply 
given two waveforms, which have arbitrary time and orientation 
offsets. They must be aligned more actively. 

We can separate this into two steps: first align the time, then 
align the frames. To align the time, we will need some measure 
of the waveform that is independent of orientation. For example, 
we can use the magnitude of the angular velocity \oj\. We then 
choose some fiducial time t&i and find the value of St such that 

\to A (ted)\ = \b> B (t M + $t)\ . (21) 

The time coordinates of waveform B may then be shifted 
as t i-» t + St. The main limitation with this method is 
that the magnitude |a)| is not always strictly monotonic for 
highly precessing systems. Usually it is possible to find a time 
for which it is monotonic. Alternatively, we can find St by 
minimizing the squared difference between the two sides of 
Eq. (21) integrated over some significant span of time. There 
may also be quantities other than |a)| that may be slightly 
more robust against this non-monotonicity or numerical noise — 
quantities such as flux or the total power in the waveform. 

Now, once the time coordinates have been properly aligned, 
it is a simple matter to align the frames. We simply apply the 
transformation R B {t) i-> RA(fgd) Rait), using Ra as given in 
Eq. (19a). Then, at ffid, there will be precisely no difference 
between the frames; in particular, <E>A(ffid) = 0- Reference [32] 
suggested essentially this same transformation, but included 
an additional rotation about the z axis because there was still 



If needed, the orbital evolution can be further isolated from the precessional 
dynamics using the minimal-rotation frame. 



one degree of rotational freedom in that paper. Here, we 
have assumed that the orientation of each waveform has been 
completely fixed at fjd, as discussed in Sec. IV A, though we 
may require a fixed rotation to the physical system of one. In 
particular, following the discussion at the end of Sec. IV A, we 
see that this transformation actually rotates a> B by Ra- Note 
that no rotation of the waveform modes is to be done here; we 
are only changing how we think of the frame in which those 
modes are decomposed. 

An alternative approach involves simultaneously fixing the 
time and frame. In previous work with nonprecessing systems, 
this was done by minimizing the squared difference in some 
quantity (A<p 2,2 , for example) integrated over some span of time. 
This was used in an effort to nullify spurious effects such as junk 
radiation or residual eccentricity [81, 82]. A similar program 
can certainly be applied to <f> A by minimizing 

T(St, Rg) := f 2 4 llog [R A (t) R B (t + St) R s f dt . (22) 

This requires simultaneously optimizing over the time offset 
and all three degrees of rotational freedom. In particular, a 
simplification that occurs in the nonprecessing case and allows 
the problem to be reduced to one dimension [83] will not work 
in the precessing case due to noncommutativity of rotations; the 
problem must remain truly four-dimensional. Nonetheless, the 
minimization is straightforward, and the frame is adjusted as 
Rfl(0 >-* RgRB{t + St). Again, the waveform modes themselves 
are not to be rotated; just the frame information. 

Once the waveforms are aligned, it is a simple matter to 
hybridize them with a slight generalization of the standard 
method. Typically, we use only information from waveform 
A before some time t\, and only information from waveform 
B after some ?2, with a transition in between. The waveforms 
are assumed to have been aligned somewhere in the range 
between t\ and t^. The transition may be accomplished with 
some (usually smooth) monotonic function r(f) that equals 1 
before t\ and after ti- Then, the hybrid version of any mode 
jljn ma y k e defined as a simple linear interpolation between 
the two waveforms: 15 

■Crid : = /I' '"(0 r(t) + f e B "\t) [ 1 - T(t)] . (23a) 

Again, however, each waveform comes with its own frame, and 
these frames have to be hybridized. As suggested in Ref. [32], 
this can be accomplished with a form of linear interpolation 
defined for rotors: 

Rhyb„d := L(r(t); R A (t), R B (t)) , (23b) 

where the interpolant L is given by Eq. (A30). As discussed 
in Appendix A 4, it is critically important to use the correct 



5 In previous work, this formula is usually applied separately to the phase 
and modulus of the mode. With this new description, there seems to be no 
advantage in decomposing the waveform in this way. The formula given 
here is written assuming complex mode data. 
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interpolant. Finally, we must note we only have the frames R4 
and Rg sampled at discrete (essentially arbitrary) points, so 
we will need to interpolate between those points to the desired 
t. For this, smoother interpolation is required. Appendix A 4 
discusses a method using cubic splines reinterpreted for rotors. 

In the computer code included among this paper's an- 
cillary files, waveform objects may be constructed with 
GWFrames. Waveform. They may be aligned, compared, and 
hybridized with methods such as AlignTime, AlignFrame, 
AlignTimeAndFrame, Compare, and Hybridize. 

V. CONCLUSIONS 

The angular velocity of a waveform was defined in Sec. II 
by the rotation which minimizes the time dependence of the 
waveform. This fairly nebulous criterion was reformulated 
precisely, and led to a simple formula, providing us with 
a geometrically meaningful description of the motion of a 
waveform. We also saw that co can be considered to be that 
vector which makes the action of the operator -i co ■ L as equal 
as possible to the action of d t , in a sense that can also be made 
surprisingly precise. 

The angular-velocity vector and the dominant eigenvector 
Vf of (LL) proposed by O'Shaughnessy et al. [31] provide us 
with powerful tools to understand and manipulate waveforms, 
with no reference to meaningless gauge quantities. Section III 
showed that these two vectors can be used very effectively and 
accurately to find at least part of the solution to the important 
inverse problem. Determining the three remaining degrees of 
freedom is left for future work, though some suggestions were 
made for how to do this. 

Beyond this fundamental benefit, co also provides key 
practical advantages. We can readily calculate the corotating 
frame, which also has angular velocity co. Transformed to this 
frame, the waveform is literally as constant as possible. When 
functions are slowly varying, they are easily approximated by 
low-order functions; they can be numerically interpolated and 
differentiated quite accurately; and fewer data points are needed 
to record their values than for quickly varying functions. This 
type of technique has already seen great success in numerical 
simulations themselves [84, 85]. 

Putting these together, we can also perform all the standard 
manipulations needed for waveform analysis. Data collected 
from a simulation at different radii can be extrapolated nicely. 
Two waveforms (e.g., different resolutions of a numerical 
simulation, or an NR and a PN waveform) can be aligned, 
compared, and hybridized readily. The only additional steps 
necessary are comparison and hybridization of the frames, but 
these are easily achieved using formulas given by Eqs. (19a) 
and (23b). Notably, the use of quaternions vastly improves 
numerics and allows us to access the geometrically meaningful 
elements of rotations. 

The code included with this paper implements all the 
techniques discussed above, showing that they are ready to 
use in waveform analysis. There are, however, issues that 
may benefit from further investigation. As mentioned, more 



work is needed to complete the solution of the inverse problem. 
Also, it is certainly possible that different techniques could 
further simplify the ringdown, for example. While preliminary 
results show that reasonable, smooth results for co and Vf are 
obtained throughout the inspiral, merger, and ringdown, the 
waveform in any frame still has very complicated structure 
during ringdown, presumably stemming from the difference 
between spin-weighted spheroidal harmonics and spherical 
ones [86, 87]. More specific methods [29, 87, 88] will likely 
be needed to adequately capture features of general ringdowns 
with simple models. 

Nonetheless, we can conclude that co and Vf already deliver 
a complete system for waveform analysis. When implemented 
with quaternion methods, the system is robust enough to be 
applied blindly to both precessing and nonprecessing systems. 
This consistency simplifies the production and analysis of both 
types of waveform. 
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Appendix A: Quaternions and rotations 

Unit quaternions clearly constitute the representation of 
choice when computing with spatial rotations. Quaternions 
have become the dominant technique in fields as diverse as 
computer graphics, robotics, molecular dynamics, navigation, 
and orbital mechanics. They are closely related to the axis- 
angle formalism, which gives us clear geometric intuition and 
avoids the problem of gimbal lock associated with singularities 
of the Euler angles. But there is also a clear notion of 
them as operators, giving us all the advantages of the matrix 
representation of rotations. They are trivially inverted and 
easily composed, and the logarithm and exponential functions 
are easy to evaluate, presenting further advantages over all other 
representations. 

For all these reasons, this paper and the accompanying code 
use the notation of quaternions. Here, the basic elements 
of quaternion math are summarized, Wigner's T> matrices 
and the spin-weighted spherical harmonics are expressed 
directly in terms of quaternions, and formulas for the linear 
interpolants and splines of rotors are given. The computer code 
included among this paper's ancillary files contains all of the 
quaternion functions discussed here. The fundamental object 
is the GWFrames. Quaternion, which has numerous methods. 
See the documentation for more details. Also, note that 
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TABLE I. Quaternion notation 


n 


V,/ Lid LCI 


Ha 


1 nmn/inunt /t/ nf trip nnQturninn 
X^OLLLpOIlCIlL Lt Ol LI1C l^UdLCIIllOIl 


o 


^OIlJUgdLC. Vt/Oi y2? ^3^ 


IQI 


Norm: V<?q + <^ + ^ + q\ 


9 


Vector part: (qi,q2,q?,) 


9 


Magnitude of vector part: V#j + q\ + q\ 


4 


Normalized vector part: q/q 


<) 


Logarithm: logQ 


q 


Magnitude of the logarithm: q 


ang R 


Angle of a rotation: 2 r = 2 |log R| 



Mathematica [89] returns incorrect results for logarithms of 
general quaternions, and is thus not a reliable tool for most of 
the calculations in this paper. 

1. Elements of quaternion mathematics 

A quaternion is a set of four numbers, usually denoted as 

Q = (<7o,<7i,<72,<73) =qo + q- (Al) 

We summarize the notation in Table I. The quaternions form 
an algebra, meaning that the quaternions form a vector space 
(over the real numbers), as well as a group where the product 
is defined by 



PQ = (poqo-p-q) + (poq + qoP + pxq) 



(A2) 



Here, the dot product and cross product of vectors take their 
usual meanings. Note that this product is neither commutative 
nor anti-commutative in general. The conjugate of a quaternion 
is defined as 

Q := (q Q , -q u ~q 2 , -qs) =qo~q, (A3) 
and the norm of a quaternion according to 

IQI 2 = QQ = ql + q\ + q\ + q\ = q\ + q • q ■ (A4) 

A unit quaternion is simply a quaternion with unit norm. 
Because quaternion multiplication is associative, we can find 
a useful inverse of a quaternion by taking the conjugate and 
dividing by the squared norm: 



Q_ 

IQI 2 



(A5) 



In particular, the inverse of a unit quaternion is just its conjugate. 
Note, however, that while a unit quaternion has norm |R| = 1, 
its square is not 1 in general. For example, if R = u is some 
unit vector, we have R 2 = - 1 . 

Now, given any vector v, we can define the transformation 
law 



v' = RvR, 



(A6) 



where the right-hand side involves quaternion multiplication 
with v interpreted as a quaternion with scalar part vq = 0. It 



is not hard to check that if R has unit magnitude, then this 
transformation law preserves orientation, angles, and lengths — 
and is therefore a rotation. These rotations compose in the 
natural way, and we will see below that we can construct a 
unit quaternion representing any desired rotation, which means 
that the unit quaternions form a representation of the rotation 
group. 16 

Using the product law for quaternions, we can define the 
exponential of a quaternion according to the standard power 
series: 



expQ := > — 



(A7) 



n=0 



Note that, because of the non-commutativity of quaternion 
multiplication, the usual rules of exponents do not apply. 
In particular, exp[P + Q] + expP expQ unless P and Q 
commute — which happens precisely when their vector parts 
are parallel. Given some angle 9 and some unit vector u, we 
can show that the unit quaternion 17 



R = exp — ii = cos — h u sin - 
F [2 J 2 2 



(A8) 



represents a rotation through the angle 9 about the axis « (in the 
positive sense, using the right-hand rule). This illustrates the 
connection between the axis-angle and the unit-quaternion 
representations of rotation. The factor of 1/2 needed in 
the exponential is a result of the double-sided rotation law, 
Eq. (A6). 

By inspection of Eq. (A8), we see that we can also define a 
reasonable logarithm of nonzero quaternions: 18 



log Q := log |Q| + - arctan — 
q qo 



(A9) 



Note that the logarithm of a unit quaternion will be a pure 
vector — log |Q| = 0. For compactness, we define the notation 
C| := log Q and q := |q|. As with the usual complex logarithm 
and the real arctangent function, this function is multivalued; 
the magnitude of the vector part is ambiguous up to integer 
multiples of 2 n. We typically choose the principal value so 
that the norm of the vector part is in [0, n], as with the complex 
logarithm. Choosing the branch must be done carefully, in order 
to obtain correct geometric results and reasonably continuous 
functions of time. When differentiating the logarithm (as in 
Sec. A3 for example), we will treat the function as being 



16 In fact, because of the double-sided rotation law, Eq. (A6), R and -R 
represent the same rotation, so the unit quaternions provide a double cover of 
the rotation group SO(3); the group of unit quaternions is actually isomorphic 
to SU(2). Unsurprisingly, the logarithms [defined in Eq. (A9)] of unit 
quaternions form a group isomorphic to su(2). Hence the notation r = log R. 

17 More generally, the exponential of any quaternion is exp Q = 
exp |Q| exp(Q/ |Q|), where the second factor can be evaluated according 
to Eq. (A8). 

18 Again, this expression generalizes the more familiar complex relation log z = 
log \z\ + i arctan 5z/iRz, where i is replaced by a general three- vector. 
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continuous. On the other hand, sometimes in the very same 
formula, we will assume the logarithm takes on its principal 
value. 

The principal value of the quaternion logarithm can actually 
be restricted further if our purpose is only to cover SO{3). 
Because of the double-sided rotation law of Eq. (A6), the final 
vector is invariant under R i-> -R. This is equivalent to 



are required and fewer singularities are encountered. This 
expression is very robust and deals well with finite numerical 
precision. This expression is ill defined whenever u + w = — 
which is not surprising, as there are infinitely many "shortest" 
ways to rotate a vector into its opposite. These are two ways of 
expressing the fact that there are infinitely many square roots 
of -1 among the unit quaternions. 



r m> 



(A10) 



In particular, we can apply this formula when r > n/2, ensuring 
that r e [0, n/2]. The transformation gives rise to different 
rotors, but the same rotation — at least for objects of integral 
spin weight, like vectors. We will use this fact when integrating 
the angular velocity in Sec. A 3 to eliminate an edge case and 
improve numerical behavior. 

As with exponents of real numbers, we can define 



exp[P logQ] 



(All) 



This formula will be usually be applied in cases where P is a 
pure real number, though other formulas may be advantageous 
in such cases — as illustrated in the case of the square-root 
below. 

The square root of a quaternion is particularly useful in 
constructing rotations taking one vector into another as directly 
as possible. We can find a formula for it with an elegant 
geometric interpretation and important numerical advantages 
over Eq. (Al 1) with P = 1/2. The product of two unit vectors 
-u w is a rotation in the ii-w plane of twice the angle between 
those vectors, in the sense from w to u. The square root of 
this product is the same rotation through only half that angle — 
in particular, V— w w is the most direct rotation taking w into 
u . We need to bisect the angle between them, and a familiar 
geometric construction that achieves this is the diagonal of the 
rhombus having w and u as sides: 



u + w 
\u + w\ 



(A12) 



Then the rotation we want is 



; „ „ UW-l l-UW ,...„, 

V-uw = +vw = + = ± — - (A13) 

I" + w \ V 2 t! -(m«0o] 

Computing the square root using this expression is easier than 
using Eq. (All), in the sense that no transcendental functions 



2. Formulas for rotations and SWSHs 

We now express Wigner's X) matrices and the spin-weighted 
spherical harmonics (SWSHs) directly in terms of quaternions, 
so that no conversion to or from the more usual Euler-angle 
representation is necessary. In the following, we will treat the 
general case in which the spin weight s is arbitrary; the formulas 
given here do not assume s = ±2. 

The SWSHs form a basis for spin-weighted functions on the 
sphere [90-92]. Goldberg et al. [49] showed that the SWSHs 
can be expressed as special cases of Wigner's D matrices, so 
that by constructing 35,,,',,,, we will obtain s Y( t „,. Defining the 
parts of the quaternion Q as 

Q a :=qo+iq3 and Q b := q 2 + i q\ , (A14) 
we can express quaternion multiplication as 

(PQ)a = Pa Qa ~ H Qb , (A15a) 
(PQ)b = Pb Qa + PaQb- (A15b) 

Quaternions are isomorphic to (Pauli) spinors, and the two 
parts of the quaternion defined here are essentially the two 
components of the spinor. The choices of signs in Eq. (A14) 
are — to some extent — arbitrary conventions. However, care 
must be taken to ensure that the resulting 35 matrices form 
a representation of the rotation group rather than an anti- 
representation, and to ensure that the handedness of space is 
preserved. Our purpose in choosing these particular signs is to 
reproduce the standard SWSHs as special cases. In particular, 
note that the presence of in the definition of Q a is what picks 
out the z axis as the point of reference on the sphere, so that 
the polar angle is measured with respect to it, rather than the x 
or y axes. 

Now, following the standard derivation [93], we obtain 



dm' —i 



r>2m 
> R b 



(-D 



{+m 



(R) 



p 



2m 



4 



(t+m)\({-m)\ 
(t+m')\({-m')\ 



\Ra\ 



26— 2m j^m+m' pm-m 



Z P (-D P ( 



e+m 
p 



\tt-m' Wl^!\ 2 ^ 
) \C-p-m) \ \R a \ ) 



when R a - 0, 
when Rb - 0, 
otherwise. 



(A16) 



This expression is valid for all integral and half -integral values of I > 0; naturally, we only need integral values I > 2 for 
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the s — -2 fields discussed in this paper. Note in particular 
that (-R) = (-l) 2m D (e) , (R), which may be nontrivial 
when to can take half-integral values. This is another statement 
of the fact that rotation through 2n may not return fields of 
half-integral spin weight to their original values, but will return 
fields of integral spin weight. This fact will be useful below 
when we integrate the angular velocity. 

To recover the usual expressions for D in terms of Euler 
angles, we use R = e"^ 2 e 6 ^ 2 e r z ^ 2 , from which we can easily 
find 



R„ = cos - e 1 
2 



Rb = sin - e 1 
2 



(A17a) 



It must be emphasized, of course, that evaluating Eq. (A 16) 
directly is faster and deals with numerical-precision issues 
better than using the form with sines and cosines. 

Now, to express the SWSHs in terms of these £ matrices, 
we adopt conventions to agree with Ref. [94], which attempts 
to establish uniform conventions for use in numerical relativity. 
We have 

,YU* <P) = (-D' S ^ l2 ) ■ ( A18 > 

These functions are implemented in the ancillary files as 
GWFrames.WignerDMatrix and GWFrames.SWSH. 



3. Integrating the angular velocity 

In many contexts, quaternion-valued functions of time turn 
up. These may be differentiated or integrated with respect 
to time, much as vector-valued functions may be. However, 
noncommutativity leads to certain problems. In the next section, 
we will see how interpolation can be handled sensibly. Here, 
we prove a vital formula used in the main text of this paper to 
integrate the angular velocity vector to find the rotor describing 
the frame with that angular velocity. 

First, we need formula for the derivative of the inverse, which 
can be obtained by differentiating QQ 1 = 1: 



V = -Q^Q- 

df v df v 



(A19) 



This is the crucial relation that allows us to calculate the angular 
velocity m of a frame described by the rotor R(f) [32]. Suppose 
that a vector vq is stationary in the rotating frame. Then, that 
vector is given in the inertial frame as v(t) = R(f) vq R(t). We 
also know that dv/dt = tttxd. Using the definition of quaternion 
multiplication and the usual commutator (Lie product), we can 
calculate m x v = v]. Another way of writing this is 

^(Ro R)= ^[w,RdoR] (A20a) 
= [RR,Rd R], (A20b) 

where the second line comes from simply evaluating the left- 
hand side and using Eq. (A19). It is not hard to show that RR 
is a pure vector, using Eq. (A 19) and the fact that any arbitrary 



quaternion Q is a pure vector if and only if Q = -Q. Then, if 
Eq. (A20) is to be true for all vectors vq, we must have 



m = 2RR 



(A21) 



The factor of 2 appears here because we are using quaternions; 
this factor does not appear in the equivalent result for rotation 
operators. 

As explained in Sec. IV, we need an expression for the right- 
hand side in terms of logarithms. To borrow notation from the 
theory of Lie groups, we define the adjoint operator using the 
familiar commutator: 



adpQ := [P,Q] =PQ-QP 



(A22) 



This notation is convenient because we will need repeated 
applications of the commutators. For example, adpQ = 
[P, [P, Q]]. Now, if P and Q are unit quaternions, their 
logarithms will be pure vectors: log P = p and log Q = q. 
We will also use the notation p := |p|, etc. Again, we can use 
the definition of quaternion multiplication in Eq. (A2) to see 
that [p, q] = 2 p x q, which allows us to use familiar properties 
of the cross product to calculate 



ad';q = 



(A23) 



(_l)C»-i)/2 [p >q ] ( 2p)»-i 

(_l)(«-2)/2 [p>q] ] (2p) »-2 

The proof is a simple induction. A standard formula [95] says 



n = 0; 
n odd; 
n > even. 



e'qe-" 



DO 1 



while a somewhat less-standard formula [96] gives us 
de» / ,„dp 



dt 



Jo 



dt 



ds 



(A24) 



(A25) 



for p < n. We can multiply this formula on the right by e~ p , and 
substitute using Eq. (A24). We then separate the resulting sum 
into three parts, corresponding to the three cases in Eq. (A23). 
These can be readily evaluated, yielding simple trigonometric 
functions, which can then be integrated:' 9 
>1 



RR= f e !t -e !t d 5 
Jo dt 

2 n. 



X 
X'( 



ad';t 



ds 



\n=0 

* ' sin(2sr) 
x + — [r, v] + 

2r 



sin (s r) 

2r 2 



(A26a) 
(A26b) 

[r, [r, t]]j ds 

(A26c) 



r + 



sin 2 r 



2r 2 



[t, i] + 



■ sin r cos r 



4r 3 



[t,[t,t]]. (A26d) 



19 Again, note the assumption that r < n, which is essential to the correctness 
of Eq. (A26), where the actual magnitude t is used. Nonetheless, we also 
assume that the derivative r exists and is continuous everywhere, which must 
be enforced by removing branch-cut discontinuities before differentiating. 
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As discussed in Sec. A 1, we evaluate the derivative r by treating 
r as a continuous function, removing any branch cuts. On the 
other hand, when used without differentiating, we have assumed 
that r < 7T. 

We can re-express this relation as a matrix equation by 
defining 



A := 



1 





0) 





1 





10 





u 






r - 



sin 2 r 



sin r cos r 





it 2 



-*3 





+ r 2 
-ri r 3 







r 2 + r 2 

~*2 *3 



-r 2 c 3 

1 7 

rf + r 2 , 



(A27) 



in which case we have the much more compact formula 

RR = y = A r . (A28) 

The determinant of the matrix simplifies to sin 2 r/r 2 , and is thus 
invertible for r < n. For the rare edge case with exactly r = n, 
the rotation expr = -1, which corresponds to the identity 
rotation, so we should have r = trr/2. For all other cases, we 
can invert the matrix explicitly to find 20 

/ r (r • ur) \ r cot r r (r • vr) 1 
r ^_A_Jj_ + J_J + _ CTXr . (A29) 

Thus, we are left with an ordinary differential equation to 
solve for r, as discussed in Sec. IV. Finally, we obtain (up 
to the constant of integration) R = exp r. We can improve 
the numerics and avoid the edge case with r = n by using the 
mapping of Eq. (A 10) between time steps whenever r > n/2. If 
the resulting rotor function is to be interpolated (or used for any 
other purpose for which R and -R are not equivalent), it may 
be useful to go back and make R(f) as continuous as possible 
by flipping the sign of R(f,) whenever |R(r,-) - R(fj_i)| > V2, 
for example. The complete algorithm is implemented in the 
ancillary files as GWFrames.FrameFromAngularVelocity. 

4. Interpolation 

When comparing waveforms, one of the most basic re- 
quirements is the ability to interpolate. The description of 
a gravitational waveform has now expanded to include both 
the SWSH modes of the waveform and the rotor describing the 
frame of that decomposition. So we need a way to interpolate 
rotors. But interpolation of rotors is complicated by the fact 
that the interpolant needs to remain normalized to unity at all 
times. While it is possible to simply interpolate the quaternions 
in R 4 and normalize the result, the interpolant will generally 
exhibit unnatural accelerations between the interpolated points, 



even in the simplest case of uniform rotation. Interpolation 
of rotation matrices is just as bad. It goes without saying, of 
course, that interpolation of Euler angles leads to complete 
nonsense — the result is highly sensitive to the orientation of the 
coordinate basis, and depends very strongly on the conventions 
for which directions the successive Euler rotations take. A 
reasonable suggestion might be to interpolate the logarithms of 
the rotors and exponentiate the interpolant. However, this also 
leads to unnatural behaviors in fairly simple cases, whenever 
the logarithms of the rotors are not parallel. Fortunately, there 
are well-motivated solutions to the problem of quaternion 
interpolation that can give reasonable results in very general 
cases. 

Recognizing that the unit quaternions can also be regarded 
as points on the unit sphere S 3 , we might further expect an 
interpolant to follow the geodesic between two points on the 
sphere. In fact, achieving this property is actually quite simple, 
using the fact that the quaternions operate as a (Lie) group. A 
simple interpolation between unit quaternions Ro and Ri that 
preserves the normalization is given by [98] 

L(t; Ro, Ri) = (Ri R ) r Ro = Ro (Ro Ri)' ■ (A30) 

Obviously, L(0; R , Ri) = R and L(l; R , Ri) = Ri, and the 
norm of L(t; Ro, Ri) is always 1. This formula is strongly anal- 
ogous to the formula for standard linear interpolation, except 
that multiplication by t becomes exponentiation and addition 
becomes multiplication. 21 This interpolation is referred to as 
"slerp" for spherical Zinear interpolation. It will be useful to 
note that 

— L(t;Rq,Ri) = log(RiRo) L(r, Ro.Ri) , (A3 la) 



= Uj; R ,Ri) log (Ro Ri) . (A3 lb) 



This formula shows us that the speed along the slerp path is 
constant, as it must be for a geodesic. 

Many problems only call for a linear interpolation like slerp. 
In particular, when blending PN and NR waveforms, each 
waveform possesses its own frame. To transition between the 
two waveforms, we must transition between the frames, which 
is just a simple linear interpolation at each instant of time where 
the extent of the interpolation (the r argument to the function 
above) depends on the time. However, we also need to be able 
to interpolate each individual waveform as a function of time. 22 
For this, we cannot use linear interpolation for the motion of the 



20 An equivalent formula was found through a very different derivation by 
Grassia [97]. 



21 This analogy should not be carried too far because quaternion multiplication 
is noncommutative. In particular, it is crucial to note that the right-hand 
side of Eq. (A30) is not equal to Ri R ~ r , for example, whenever Ro and Ri 
do not commute; such a formula actually gives very poor interpolation in 
many cases. Equation (A30) is preferable because the path it describes is a 
geodesic in the space of unit quaternions. 

22 For example, the PN waveform and the NR waveform will generally be 
calculated at different instants of time. To compare them, we need to be able 
to interpolate the values of one waveform onto the time steps of the other. 
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frame of either waveform. If we did, we would see the frame Unfortunately, the various methods — while being equivalent 
abruptly change rotation speed as it goes through each original for real numbers — are not equivalent when using quaternions 
data point, just as a linearly interpolated graph changes slope because of noncommutativity [99]. It is not clear that any 



abruptly as it passes through each original data point. Instead, 
we would prefer some higher-order technique. 

We can approach this problem in analogy with the construc- 



particular formulation will give better results than any other, so 
we may take the pragmatic approach of simply choosing one 
which is easily implemented. The result will be a spherical 



tion of curves in space, which suggests various approaches such interpolation based on the quadrilateral of a standard spline, 
as the de Casteljau algorithm for constructing Bezier curves. referred to as squad. 



In that spirit, we will define the cubic-spline interpolant in terms of the linear interpolant: 

C(f;R/,A ; ,B,- +1 ,R /+1 ) = l{2ti{\ — t,); Uji, R,-, R;+i), Uji, A;, B !+ i)) , 



(A32) 



where A, and B, + i are "control points" to be solved for. We have also defined r,(f) = (f - f,)/(f,+i - f,), where i is assumed to be 
the index of the nearest time sample such that f, < t. We can evaluate the derivative 

(A33a) 
(A33b) 



— C = — {exp[2r/(l - t,) log (l(t;; Aj, B i+ i) L(t ,; R;, R,+i) -1 )] L(t;; R„ R,+i)| 

= (2-4t,-) log(L(T,;A ! -,B, +1 )L(r i ;R,,R, +I r 1 )c + 2r,(l -t,)G + C log(R,R, +1 ) 



where G is a complicated expression, which is easy to compute, but messy to write; fortunately do not need to evaluate it because 
that term drops out when we evaluate at r* = or t,- = 1, as the factor in front of G goes to zero. We wish to ensure that the time 
derivatives are equal at the end of one segment and the beginning of the next: 



dr 



C(r,_ 1 ,R|_i,A j _ 1 ,Bi,Bi) 



= — C(t,-, Rj, A^Bf+i, R;+i) 



T;=0 



(A34) 



Note that we differentiate with respect to f, rather than t,, to account for differences in the time steps of the given data.Plugging in 
the result of Eq. (A33) and simplifying, we get 



Ati 



\- {-2 logfBiRijRi+R; log(Ri_! R,)} = ^ {2 log (A, R,) R, + R,- log(RiR i+1 )} , 



or equivalently 



l — R, {-2 log (R,- B,. ) + log(R,_, R,)j = ^ R, {2 log (R,- A,) + log(R ; R, +1 )} 



(A35) 



(A36) 



We need one more condition to solve for both variables A,- and B,-. We may choose 2 ' to set the velocity at either side equal to the 
average velocity of linear interpolations on those two sides, giving us the following two equations: 



logCRj R,- +1 )/Ar,- + log(RM R,)/Af,_ 1 1 , ,o v<[ 

R; ^ = aTT ' I g ' '' '' / log ^ R '- 1 R!) ) 

log(R,- R, + i)/Af ; + log(R,-i R,)/Af,_ 1 1 , / - . \ , , s _ ,i 
Ri = — R ; [2 log (Ri A,j + log(R,- R, + i)J . 



We can now solve for the control points: 

A,- = R,- exp 
B ; = R, exp 



log(R,-R i+1 ) + log(fi,_, Ri)M7A^-i - 2 log(R, R, +I ) 



log(RiR !+ i) Ak-i/Ati + log(R/_i R) - 2 log(R,_i Rj) 



(A37a) 
(A37b) 

(A38a) 
(A38b) 



To apply these formulas to the edge cases of i — and i — N - 1, we also define the quantities R_i = RqRi Rq and R^ 
Rjv-i Rn-2 Rjv-1) which roughly represent straight-line motion. 

I 



This choice has the nice property of agreeing with our intuition in the case ^ of « straight . line " motion . To be precise: if the transformation from R,_, to 
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In the computer code included among this paper's ancillary 
files, the functions GWFrames . Slerp and GWFrames . Squad im- 
plement linear and cubic interpolations of rotors. 

Appendix B: Rotating spin- weighted functions 

Gravitational radiation is a complex field of nonzero spin 
weight, meaning that it picks up a position-dependent phase 
under rotation [92] . The reason for this is its definition with 
respect to a dyadic which is itself defined in terms of a coor- 
dinate basis; when the coordinates rotate, the dyadic rotates. 
Depending on details of the definition of the gravitational- 
wave field, the spin weight may be s — 2 or s — -2 — the 
most common choice being the latter. Throughout the rest of 
this paper, we have assumed s = -2; in order to discuss the 
properties of general spin-weighted fields, this appendix will 
apply to general values of s. 

Suppose we have a field / of spin weight s on the sphere. To 
measure this field, we first need some standard basis for our 
space, (x, y, z). We can define the usual spherical coordinates 
relative to this basis, and write the field as a function of the 
coordinates, so that in some particular direction h, we have 
/(«) = /(#, <p). We define the rotor 

R (lM :=e** /2 e** /2 (Bl) 

and note that 

n = R(,» # ) z R(#,</>) • (B2) 

Now, suppose we rotate the physical system by some R p h ys - 
Then, the corresponding direction in the rotated field is n = 
Rphys h Rphys- We know that there must be some angles Ip) 
such that 

But these conditions are not enough to fully restrict our 
rotations. There is some other angle 24 y such that 

R (£,v>) = R phys e7Z/2 ■ (B4) 

The term involving y represents an initial rotation through 
that angle in the positive sense about the direction z, which is 



R, is written as multiplication by R, R;_i, then "straight-line" motion occurs 
when Rj+i = R, R; i R, (and for simplicity, we assume Af,_i = Atf). Then 
this average velocity is R, log(R,_i R;)/Af, which is precisely the velocity 
of a linear interpolation at that point. 

This angle is required to account for the full three-dimensional freedom in 
choosing R p hy S - It is always possible to find such an angle. However, this 
angle need not be unique for certain orientations; y may be degenerate with ip. 
Similarly, because of the familiar singularities of the spherical coordinates, 
there may not be a unique choice of (??, <p) or tp) for certain positions. 
Nonetheless, the rotations R(0 # ) e rz ' 2 and R^ ~ generated by these angles 
will be uniquely determined, much as the North and South Poles are uniquely 
determined despite the ill-defined longitude at those points. 



equivalent to & final rotation about the direction n . For spin- 
weighted functions, this corresponds [92] to multiplication of 
the function value by e~ lsy . Thus, in this basis, we measure 
a different field /, related to the field / measured in the first 
basis by 

m<p) = m<p) e - is T . (B5) 

For s = 0, we recover the familiar result that a scalar field does 
not depend on the frame in which it is measured. 

The spin-weighted spherical harmonics (SWSHs) form a 
basis for spin-weighted functions on the sphere [90-92], just 
as standard spherical harmonics form a basis for spin-zero 
functions. We can write 

/(#, V ) = 2/''\y f , m (^), (B6a) 

t,m 

M 0) = Yj f f ' m ' Y i»@> ® ■ ( B6b ) 

t,m 

The SWSHs themselves are just special cases of the Wigner D 
matrices (see Eq. (A 18) and Ref. [49]). We can then use the 
fact that the D matrices form a representation of the rotation 
group to find the transformation law for SWSHs: 

^' m (Ri R 2) = 2 I>LV R i) ^-, m ( R 2) ( B7 > 

m" 

implies, using Eq. (B4), that 

t Y t J$, 0) = £ s Yw(&, <p) Xj^m<( R phys) e - is ? . (B8) 

m' 

The dependence of y on <p) means that, strictly speaking, 
the SWSHs with s + do not transform among themselves 
under rotations. Naturally, when coupled to the appropriate 
spin- weighted tensors, the complete object transforms as ex- 
pected [16]. Similarly, the modes f e,m transform nicely thanks 
to a convenient cancellation. Inserting Eqs. (B6) and (B8) into 
Eq. (B5), we can show that 

f' m = ^f Xm '^J^)> (B9a) 

m' 

or equivalently 

m' 

These are precisely the same as the transformation laws for 
modes of standard (s = 0) spherical harmonics, and do not 
depend on y. 

It is worth pointing out that a rotation of the physical system 
by Rphys is equivalent to a rotation of the basis with respect to 
which that system is measured by Rf ram e = Rphys- But there is an 
important subtlety to be observed in the context of composing 
rotations. Physical rotations compose by left multiplication, 
whereas rotations of the frame compose by right multiplication. 
That is, if we first perform a physical rotation R p i then rotate 
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that system by R p 2, it is equivalent to rotating the original 
system by R p 2 R p i — just as with vectors. On the other hand, if 
we first rotate the frame by Rn, then by Rf2, it is equivalent to 
rotating the frame by R tl Rq — which is opposite to the usual 
behavior. We must carefully bear in mind the type of rotation 
we are performing. 

Expressions for the Wigner D matrices are given directly 
in terms of the rotor in Eq. (A 16), which avoids the need 
for conversion to Euler angles. The SWSHs are expressed 
as particular components of these matrices in Eq. (A 18). In 
the computer code included among this paper's ancillary 
files, the Wigner X) matrices and SWSHs are implemented 
as GWFrames.WignerDMatrix and GWFrames.SWSH. Wave- 
form objects may be constructed with GWFrames. Waveform, 
and transformed to different frames using the methods 
RotatePhysicalSystem and RotateDecompositionBasis. 



1. Maximization 

In general, we can describe the process of finding the 



radiation axis as a maximization over R of the quantity 

,2 



Appendix C: Other methods of choosing a frame 

In the interests of completeness, and to facilitate direct 
comparisons using common language, we now review three 
other methods of choosing a frame to eliminate mode-mixing 
in waveforms from precessing systems. Each of these methods 
constructs a new frame by ensuring that the z direction lies 
along some chosen axis which is roughly the axis of rotation 
of the waveform. These differ from the corotating frame 
introduced in the main text of this paper, in that the waveform is 
still rotating in these new frames. The considerations of Sec. Ill 
(and in particular the left panel of Fig. 1 ) suggest that among 
these three, the preferred method is that of O'Shaughnessy et al. 
supplemented with the minimal-rotation condition [32]. In 
particular, when setting the integration constant discussed near 
the end of Sec. IV A, that is the method of choice. However, 
over all, the corotating frame is generally still a preferable 
choice. 

We first describe two methods suggested by Schmidt 
et al. [30] and O'Shaughnessy et al. [31], using a common 
notation which allows a common implementation by means of 
explicit maximization of a quality function. We then describe 
the method of O'Shaughnessy et al. in the way in which it 
was introduced, which allows a second implementation by 
solution of an eigensystem. A third possible axis suggests 
itself given the results of this paper: the angular-velocity vector 
a> given by Eq. (7). All three need an additional step to 
remove sharp features in the waveforms, given by the minimal- 
rotation condition. In this section, we review each alternative 
in the language of quaternions, suggesting improvements for 
numerical accuracy and robustness. 



t,m 



(Cla) 
(Clb) 



Here, the wt, m are simply weighting factors. Schmidt et al. 
took these factors to be W2,±2 = 1, and zero otherwise; 
O'Shaughnessy et al. effectively chose W( m = m 2 , with some 
cutoff { above which w^ m = 0. 

This function is actually degenerate with respect to initial 
rotations about the z axis, because such rotations simply affect 
the overall phase of the term inside the absolute value. For 
numerical efficiency, we need to restrict Q to some nondegen- 
erate domain for efficient numerical maximization. Whereas 
previous references [30-32] used rotations of the form R(# ; ^), 
we choose instead to use rotations of the form R(#^) e~^ 2 . All 
such rotations can be written as R„ = e" for some vector v in 
the x-y plane, of magnitude less than or equal to n/2. Using 
rotations of this form significantly simplifies calculation of 
the Wigner £ matrices and eliminates the degeneracy near the 
identity, which substantially improves numerical accuracy and 
stability for mildly precessing systems. 

Of course, the sphere cannot be covered homeomorphically 
by a single coordinate chart, so an additional degeneracy 
remains: all vectors on the boundary of our domain result 
in rotations with equal values of Q. However, this set has 
measure zero in the domain itself, meaning that it is almost 
never encountered. Moreover, the effect of this degeneracy 
will be completely eliminated in Sec. C4. In fact, we find it 
convenient to extend the domain further. We maximize Q(e v ) 
for all vectors v in the entire x-y plane, parameterizing the 
function arguments by the usual coordinates (x, y) eR 2 . There 
are now degeneracies on circles of radius n n/2 centered on the 
origin for all integers n > 0. Again, however, these degeneracies 
cause no practical difficulties. 

Given values for the modes f e,m and the weights W( jn , the 
right-hand side of Eq. (Clb) is known analytically, using 
Eq. (A 16), as are its derivatives with respect to x and y. These 
functions are ungainly, but can be written down explicitly, 
plugged into a computer, and used in efficient numerical 
optimization routines. Direct maximization of Eq. (CI) is 
simple to implement, and can be made reasonably efficient 
and robust. It does have disadvantages, however. At each step 
in the minimization routine, all relevant X> matrices need to 
be recomputed. When the W( m are nonzero for many values, 
this can become very expensive. In such cases, it can be 
significantly more efficient to find a radiation axis using the 
following method. 

In the computer code included among this paper's an- 
cillary files, a waveform object may be constructed with 
GWFrames .Waveform. The axis suggested by Schmidt et al. may 
then be found by applying the SchmidtEtAl Vector method. 
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2. Dominant principal axis 

In general, if w^ m = W{ m 2 where Wt only depends on I, 
then this can be presented in a different form and solved as an 
eigenvector problem — which is the approach O'Shaughnessy 
et al. actually used when introducing their method. Define 25 

(L ia L b) ) := J] w * < { > m '\ka L b) \i, m) f Cjn , (C2) 

£,m,m' 

where L a is the usual angular-momentum operator. The 
radiation axis is chosen to be the dominant principal axis 
Vf of this tensor — the eigenvector with the eigenvalue of 
largest magnitude, which can be found with standard algebraic 
techniques. We can find some rotation R ax taking the z axis 
into the dominant principal axis. Reference [32] showed that 
such a rotation maximizes the function Q of Eq. (CI). 

Again, however, this rotation is not unique. Moreover, 
the dominant principal axis is only defined up to a sign, 
and numerical implementations may choose between the two 
options effectively randomly. A naive choice of R ax (f), then, 
may flip back and forth discontinuously. Fortunately, we can 
overcome this problem easily by taking a, i-> whenever 
a, • < 0. Then, we can ensure that the appropriate R ax (f,) 
is as close 26 as possible to R ax (f,_i) by choosing 



Ra ■- y-&i fij-i , 

Rax(f/) = RARax(fi-l) 



(C3a) 
(C3b) 



In Eq. (C3a), the vectors are multiplied as quaternions, and the 
square root may be found with the help of Eq. (A 13). We start 
out with a i = z, and build up the frame by stepping forward 
in time according to Eq. (C3), where each a, is the dominant 
principal axis at that instant of time. However, for reasons 
of numerical stability, fi,_i in Eq. (C3a) should be expressed 
as R ax (f;_i) z R ax (f,_i), rather than as the principal axis at the 
previous instant. 

The frame found by this method has certain advantages over 
maximization of Eq. (CI). In particular, the matrix in Eq. (C2) 
need only be computed once for each time step. The dominant 
principal axis is then obtained from this. No calculations of 
Wigner's D matrices are necessary, which tends to make the 
computation fast. Also, some minor care is needed to make the 
present method robust, mostly involving choosing the direction 
of the axis to be consistent from moment to moment. 

In the computer code included among this paper's an- 
cillary files, a waveform object may be constructed with 
GWFrames. Waveform. The dominant principal axis of {L^Lb)) 
may then be found by applying the OShaughnessyEtAlVector 
method. 



25 O'Shaughnessy et al. used wz = l, and for all other weights, as well as an 
overall normalization which is ignored here for simplicity. 

26 The distance between two rotations Ri and R2 can be defined as 
2 |log(R[ R2) , which is the minimum angle needed to rotate one into the 
other. See Appendix A 1 for more details. 



3. Aligned with the angular velocity 

A very similar frame can be defined, using the angular- 
velocity vector a> in place of the dominant principal axis of 
(L( a Lb)). The vector co was found in Sec. II A and is given 
explicitly by Eq. (7). This can be used for the «, in Eq. (C3). 
Note that using the only the direction of the vector to align 
the axis of the new frame throws away some information. 
Specifically, the magnitude of a> is meaningful and is used 
in Sec. IV to derive a fully corotating frame. Nonetheless, 
the waveform rotation in this frame is about the z' axis at 
each instant, making the time dependence of the waveform 
in this frame quite similar to that of a nonprecessing system in 
a stationary frame. 

In the computer code included among this paper's an- 
cillary files, a waveform object may be constructed with 
GWFrames. Wave form. The angular velocity may then be found 
using the AngularVelocityVector method on such an object. 

4. The minimal-rotation condition 

Each of the three methods discussed above is critically flawed 
when applied to a time-series of data, unless followed by the 
procedure described here. The end result of any of the three 
previous methods is some rotation R ax (f) that takes the z axis 
into the chosen radiation axis: a(t) = R ax (?) zR ax (f)- As 
mentioned, however, this is by no means the only such rotation. 
Indeed, because of the invariance of z under rotations about 
the z axis, any rotation of the form 



R(f) = R ax (?) exp 



r(0 - 



(C4) 



will do the same. Arbitrarily setting y{f) = leaves us with 
large extraneous features in the phase of each mode of the 
waveform. Reference [32] showed that it is easy to impose 
a condition on y(t) such that the total rotation R satisfies a 
geometrically and physically meaningful criterion referred to as 
the minimal-rotation condition. This section simply reiterates 
the previous description in quaternion form and suggests a more 
accurate way of finding R ax in some cases. 

To motivate this condition, we first define the radiation 
frame's instantaneous angular-velocity vector -nr. Then, for any 
vector v that is stationary in the radiation frame, its derivative 
in an inertial frame is given by 



v = mxv , 



(C5) 



where a dot denotes differentiation with respect to time. A 
radiation frame is — by definition — a frame in which the radia- 
tion axis a is stationary. So, its derivative in an inertial frame 
is a = rn x a. Taking the cross product of both sides of this 
equation by a, using the standard vector triple product formula 
with the fact that a has unit magnitude, then rearranging, we 
find 



vj = axa + (a-vj)a. 



(C6) 



Now, we might hope that since a( f) is actually measured from 
the waveform, this might be enough to specify the frame. 
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Unfortunately, Eq. (C6) defines the component of w along a 
circularly; it is undetermined, so we need another condition. Of 
course, an obvious solution presents itself. When the radiation 
axis is stationary, we can expect that the frame should be 
stationary. To achieve this, Eq. (C6) shows that we must have 
tn ■ a — 0. Because a is a geometric object, independent of 
the frame in which it is measured, this relation is geometrically 
meaningful. We therefore require this condition even in the 
nonprecessing case. This minimizes the magnitude of rn, so 
we refer to it as the minimal-rotation condition [18, 32]. We 
adopt this condition as the criterion for selecting the radiation 
frame. 

Of course, the frame is not given by its instantaneous rotation 
vector, but by its orientation at each instant of time. So we need 
to express m in terms of R(f) to impose our condition. This 
is conveniently calculated in Sec. A 3, which shows that tn - 
2 R R. Because the radiation axis is given by a — R z R, the 
minimal-rotation condition becomes trr-a = 2RRRzR = 0. 
Invariance of the dot product under rotation shows that we can 
also write this as RR - z — 0. Expanding R as given in Eq. (C4), 
we see that the minimal-rotation condition is satisfied if y{t) 
satisfies 

y{t) = -2 R ax (f) R ax (f) -| = 2 (R ax (f) R ax (f) z\ , (C7) 
where the subscript here denotes the scalar part. Now, 



since R ax is assumed to be known — perhaps by one of the 
three foregoing methods — we can evaluate the right-hand side, 
then integrate in time, and insert the result into Eq. (C4). 
Note that the integration constant y(0) is undetermined. This 
corresponds to the usual freedom in choosing a phase, familiar 
from nonprecessing systems, and will have to be fixed in a 
similar way. 

As a practical matter, the rotor R ax is typically computed 
using Eq. (A 13) with w = z and u — a. We can easily 
differentiate this, assuming w is constant, and arrive at an 
analytical formula for R ax in terms of a. In the construction 
of PN waveforms, the latter is known analytically, and may be 
inserted into this formula for higher accuracy: 

R ax = ±d, l 1 ~ UZ = (C8a) 
V2[l-(az) ] 

/ -d,az d,a 3 \ 

=± (ww"^iH- (C8b) 

Here, we have used -(a z)o = «3 for simplicity. 

In the computer code included among this paper's ancil- 
lary files, an array of quaternions can be put into minimal- 
rotation form using the GWFrames . MinimalRotation function. 
A waveform object constructed with GWFrames .Waveform can 
be transformed into the frames discussed in this section using 
methods beginning with TransformTo. 
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